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ABSTRACT 


A  comparison  of  standard  z  transform  and  algebraic 
substitution  synthesis  methods  for  digital  filters  are 
presented  with  emphasis  on  the  frequency  domain  character- 
istics.  A  generalization  of  Martin's  procedure,  which  leads 
to  recursive  filters,  is  obtained.   Also  a  new  method  to 
synthesize  a  digital  filter  from  a  magnitude  squared  versus 
frequency  characteristic  specification  is  presented. 
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I.   INTRODUCTION 

Digital  filters  have  become  increasingly  important 
because  of  technological  developments  in  integrated  circuits 
and  digital  techniques.   This  progress  has  increased  the 
number  of  applications  in  which  digital  filtering  techniques 
can  be  used.   In  spite  of  the  fact  that  it  is  a  relatively 
new  development,  a  vast  body  of  literature  is  available  for 
digital  filters.   The  relationship  between  digital  filters 
and  classical  continuous  filter  theory  can  be  depicted  as 
in  Fig.  1.1. 

Start  from  the  differential  equation  at  the  center  of 
the  diagram.   The  Laplace  transform  of  the  differential 
equation  may  be  taken  and  a  transfer  function  formed.   The 
frequency  spectrum  of  the  continuous  filter  may  be  evaluated 
as  shown  in  the  diagram  by  letting  s  =  joo.   The  upper  part 
of  the  diagram  shows  the  relationship  between  the  transfer 
function  of  the  continuous  filter  in  the  frequency  domain, 
and  the  frequency  domain  characteristics  of  the  sampled  impulse 
response  of  the  continuous  filter.   In  deriving  these  results 
use  is  made  of  the  Z  transform.   The  lower  part  of  the  diagram 
shows  numerical  integration  of  the  differential  equations. 
Adopting  a  numerical  integration  algorithm  yields  a  difference 
equation  which  can  generally  be  expressed  as  a  discrete  time 
equation.   For  a  given  integration  scheme,  as  illustrated  in 
the  diagram,  it  is  possible  to  substitute  a  function  of  z 
(depending  upon  the  type  of  integration)  for  the  variable  s 
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of  the  Laplace  transformed  continuous  equation.   When  z  = 
exp  (joot)  is  substituted  in  this  equation  a  spectrum  is  obtained 
which,  in  general,  is  not  identical  to  the  spectrum  of  the 
sampled  impulse  response  of  the  continuous  filter.   One  of  the 
objectives  of  this  thesis  is  to  compare  numerical  integration 
with  the  standard  Z  transform  results  as  depicted  in  the  upper 
portion  of  the  diagram. 

The  Z  transform  technique  is  described  in  the  work  of  J. 
Ragazzini  and  L.  H.  Zadah  [1]  and  has  been  discussed  by  many 
other  authors  including  B.  C.  Kuo  [2].   J.  Salzer  [3]  has 
written  a  classic  paper  on  the  comparison  of  digital  integra- 
tors in  the  frequency  domain  with  the  ideal  continuous 

2  z-1 
integrator.   The  bilinear  substitution,  s  =  —  — r=r   /  which 
3  T  z+1 

corresponds  to  trapezoidal  integration,  was  considered  very 
early  by  A.  Tustin  [4]  and  has  been  discussed  thoroughly  by 
K.  Steiglitz  [5] . 

In  this  thesis  several  areas  of  digital  filter  theory  are 
considered  in  detail. 

Chapter  two  contains  a  derivation  of  basic  Z  transform 
theory,  a  review  of  the  synthesis  of  digital  filters  using 
standard  Z  transforms,  and  a  discussion  of  the  effect  of 
sampling  frequency  on  the  frequency  spectrum.   The  relationship 
between  the  s  and  z  plane  representation  of  the  equations  is 
presented.   Chapter  three  contains  a  general  interpretation  of 
the  algebraic  substitution,  s  =  f(z),  for  several  numerical 
integration  techniques,  a  discussion  of  the  bilinear  substitu- 
tion (trapezoidal  integration) ,  and  a  development  of  the  effects 
of  sampling  time  on  the  precision  of  numerical  calculations. 


Chapter  four  discusses  the  synthesis  of  digital  filters  from 
frequency  spectrum  characteristics.   It  contains  a  derivation, 
discussion.,  and  a  practical  application  of  a  direct  method 
of  synthesis  originated  by  Martin  [6] .   A  new  generalization 
of  this  method  to  obtain  recursive  filters  from  a  specification 
of  magnitude  versus  frequency  characteristics,  is  presented. 
The  advantages  of  this  technique  include  the  fact  that  a  zero 
phase  shift  filter  can  be  obtained.   Derivation  of  another  new 
method  of  finding  the  coefficients  of  a  digital  filter,  starting 
from  a  magnitude  squared  criterion,  is  presented.   The  advan- 
tages of  this  technique  include  reducing  computer  storage 
requirements,  and  reducing  the  number  of  multiplications  to  one 
half  of  those  required  by  the  Martin  method.   Chapter  five 
presents  a  summary  and  several  questions  for  future  research. 


II.   THE  Z  TRANSFORM  TECHNIQUE  APPLIED 
TO  DISCRETE  TIME  SIGNALS  AND  DIGITAL  FILTERS 

Systems  may  be  classified  as  linear,  time  invariant, 
causal  and  so  on.   They  may  also  be  classified  according  to 
the  nature  of  the  signals  present.   A  continuous  signal  is 
expressed  as  a  continuous  (or  piecewise  continuous)  function 
of  the  independent  variable,  t.   The  signal  must  be  uniquely 
defined  at  all  values  of  t  within  a  given  range,  except 
possibly  at  a  denumerable  set  of  points,  as  for  example  in  a. 
square  wave.   A  discrete  signal  is  defined  only  at  a  sequence 
of  discrete  values  of  the  independent  variable  t.   A  quantized 
signal  is  one  which  can  assume  only  a  denumerable  number  of 
different  values,  but  quantized  signals  may  be  either  continuous 
or  discrete,  as  discussed  by  P.  M.  DeRusso  et  al,  [7]. 

There  are  many  examples  of  continuous  systems  in  nature, 
and  they  are  to  be  found  almost  everywhere.   Examples  of 
discrete  and  quantized  systems  are  perhaps  more  difficult  to 
conceive  but  this  does  not  mean  they  are  scarce..   For  instance, 
the  output  of  a  continuous  system  when  sampled  for  digital 
computation,  or  other  purposes,  becomes  a  discrete  system. 
Some  systems  are  discrete  by  nature.   The  distance  to  an  object 
as  measured  with  a  pulsed  radar;  or  a  pulse  coded  communica- 
tions system  which  is  time  multiplexed  are  examples  of  inherently 
discrete  physical  systems.   A  digital  computer  is  an  excellent 
example  of  a  discrete,  quantized  system.   The  quantizing  levels 
are  determined  ultimately  by  the  finite  number  of  bits  of  the 


number  representation.   The  discrete  characteristic  is  deter- 
mined by  the  fact  that  the  independent  variable  assumes  values 
at  finite  increments  only. 

In  order  to  deal  with  discrete  time  quantized  systems, 
existing  Z  transform  theory  has  been  extended  [Ref »  2]  ..   This 
theory  was  developed  primarily  for  sampled  data  control 
systems  and  resembles  Laplace  transform  theory  for  continuous 
systems  in  that  it  is  operational  in  nature.   The  Z  transform 
has  proven  to  be  useful  for  the  study  of  all  discrete  time 
systems  and  it  is  extended  approach  which  is  reviewed  in  this 
chapter.   Several  important  fundamentals  are  developed  and  a 
basic  misapplication  of  the  theory  which  occurs  commonly  in 
the  literature,  is  discussed. 

A.    DEFINITION  OF  THE  Z  TRANSFORM  OF  A  DISCRETE  SIGNAL 

Consider  a  continuous  signal  x(t),  for  t>0,  which  is 

sampled  at  a  frequency,  f   or  equivalently ,  every  T  =  I/f 

s  s 

seconds . 

Assume  that  the  signal  is  sampled  ideally,  that  is~,  the 
sampled  signal  assumes  the  value  of  the  continuous  signal  with 
infinite  accuracy  at  discrete  values  of  the  independent 
variable,  t  =  nT;  n  =  0,1,2,  ... 

The  process  of  sampling  can  be  described  as 

x*(t)  =  x(t)  s(t)  (2.1) 

where 

x*(t)  =  the  sampled  signal 
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x(t)  =  the  continuous  signal 

s (t)  =  the  sampling  function,  representing  a  train  of 
impulses,  T  seconds  apart. 

00 

s(t)  =  1  6  (t-nT)  (2. 2) 

n=0 

From  (2.1)  and  (2.2) 


x*(t)  =  Z    x(nT)  6  (t-nT)  (2.3) 

n=0 


Taking  the  Laplace  transform  of  the  sampled  signal  we 
have  from  Eq.  (2.3) 

X*(s)  =  £  (x*(t) } 

00 

=  Z    x(nT)  £ {6 (t-nT) }  (2.4) 

n=0 

But   £{6  (t-nT)}  =  e"nTs  (2.5) 

so  that  Eq.  (2.4)  becomes 

00 

X*(s)  =  Z    x(nT)  e"nTs  (2.6) 

n=0 

Now,  define  a  new  variable  z,  such  that 

z  =  esT  (2.7) 

so  that  Eq.  (2.6)  becomes 

oo 

X(z)  =  X*(s)        =  S    x(nT)  z~n  (2.8) 

sT    n=0 
z=e 

The  standard  procedure  to  find  the  Z  transform  of  a 
discrete  signal  can  be  summarized  as  follows: 

11 


(a)  Sample  the  continuous  signal 

(b)  Find  the  Laplace  transform  of  the  sampled  function 

sT 

(c)  Substitute  e    =  z  in  the  Laplace  transform  of  the 

sampled  signal  to  obtain  the  Z  transform. 
Using  common  notation,  the  above  can  be  written  as  follows: 


X(z)  =  £  {x(t)>  =  [£{x*(t)}] 


(2.9) 


z=e 


sT 


This  process  of  obtaining  the  Z  transform  of  a  signal  is 
illustrated  in  Fig.  2.1. 

1 .    Obtaining  the  Z  Transform  in  Closed  Form 

The  foregoing  definition  of  the  Z  transform  provides 
a  method  for  obtaining  the  transform  of  a  signal,  but  it  has 
two  disadvantages: 

a.  Results  are  presented  in  open  form  because  the 
transform  of  a  signal  is  given  as  a  summation  of  inverse 
powers  of  z  whose  coefficients  are  the  values  of  the  signal  at 
the  sampling  times.   For  sampled  functions  derived  from  contin- 
uous functions  with  poles  in  the  left  half  s  plane,   such 
coefficients  will  eventually  decay  to  zero.   If  the  decaying 
process  is  slow  compared  to  the  sampling  interval,  a  large 
number  of  terms  in  the  sampled  Z  transform  expansion  are  required 
to  represent  the  successive  values. 

b.  The  process  of  obtaining  the  Z  transform  may 
be  very  long,  particularly  when  only  the  Laplace  transform  of 
the  continuous  signal  is  known.   The  process  requires  taking 
the  inverse  Laplace  transform  of  the  continuous  signal,  sampling 
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the  resulting  time  function,  and  finally  expressing  the  sample 
values  as  a  series  of  inverse  powers  of  z. 

There  is  an  alternate  method  for  obtaining  the  Z 

transform  of  a  signal  that  is  particularly  suitable  when  the 

Laplace  transform  of  the  continuous  signal  is  known.   This 

latter  approach  gives  the  result  in  closed  form.   From  Eq. 
(2.1) 

X*(s)  =   £  (x(t)  s(t) }  (2.10) 

and  so   X*  (s)  =  X (s)  *  S(s)  (2.11) 

where  the  asterisk  denotes  convolution. 


c+j°° 

X(A)  S(s-A)  dX  (2.12) 

c-j°° 


=  2¥J  j 


now  S ( s ) 


f  00   oo 

E   6  (t-nT)  e~St  dt 
o  n=0 


00 

o  /  \     v     -nTs 
S (s)   =  Z    e 

n=0 


S(s)   =  1  /  (l-exp(-sT))  (2.13) 

provided  |exp(-nT)|<  1 

In  (2.11),  use  has  been  made  of  the  sifting  integral 
property  of  the  impulse,  namely 

6  (t-nT)  f(t)  dt  =  f (nT)  (2.14) 

'  o 
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Therefore,  the  Laplace  transform  of  the  sampled 
function  can  be  written  in  terms  of  the  complex  convolution 
integral 


X*(s)    = 


2tt  j 


fC+j00 

X(X) 
c-j00 


1-e 


-(s-X)T 


dX 


(2.15) 


If  T.  as  shown  in  Fig.  2.2  is  taken  as  the  path 
enclosing  the  left  half  of  the  s  plane  then  integral  (2.15) 
may  be  divided  into  two  integrals  cs  follows 


X*(s)  = 


V        X(X) 

1 


2-iTJ  ;  r,       1-e 


-(s-X) 


dX 


2irj 


C  X(X) 


1-  e 


(s-X) 


dX 


(2.15) 


Note  that  r..  is  composed  of  the  straight  line  joining 
the  points  (c,-j°°)  and  (c,+j°°)  and  the  counter-clockwise  semi- 
circle or  infinite  radius. 

If  the  condition 


lim   s  X(s)  =  0 

s  ■*■   °° 


(2.17) 


is  met,  then  the  second  integral  tends  towards  zero  and  Eq. 
(2.16)  becomes 


X*(s)  = 


2ttj 


X(X) 


1-  e 


-(s-X) 


dX 


(2.18) 


A  very  important  case  occurs  when  the  Laplace 
transform  X(s)  of  the  continuous  signal  is  a  ratio  of   two 
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two  polynomials 

X(s)  =  §{§}  (2.19) 

In  this  case  the  function,  X(s),  has  a  finite  number 
of  poles  and  the  integral  around  the  closed  path  can  be 
evaluated  using  the  residue  theorem.   Thus  Eq.  (2.18)  becomes 


X*(S)  =  E    ^ni 1  (2.20) 

n=l  D' (X    )  1-  e  ls  AnJi 
n 

where  k  is  the  number  of  simple  poles  of  X(A)  located  at  X    , 

with 


=   d  D(8) 

n       ds 


(2.21) 


n 

and  substituting  Eq.  (2.7)  into  (2.20),  yields  X(z)  the  Z 
transform  of  x(t) 

k    N(A  ) 

X(Z)  =  I    ^ ±-= ('2.22) 

n=l  D' (X    )  1-  enz_1 
n 

The  above  formula  (2.22)  represents  a  closed  form 
for  the  Z  transform  of  the  signal  which  contains  as  many  terms 
as  the  degree  of  the  denominator  polynomial  of  the  function 
X(s).   Expansion  of  Eq.  (2.22)  in  a  series  of  inverse  powers 
of  z  yields  Eq.  (2.8). 

B.    SYNTHESIZING  A  DIGITAL  FILTER  WITH  THE  STANDARD  Z  TRANSFORM 

In  the  previous  section  two  methods  have  been  developed 
to  obtain  the  Z  transform  of  a  sampled  signal,  starting  from 
either  the  continuous  time  function  or  its  Laplace  transform. 

17 


The  same  procedure  can  be  followed  to  obtain  the  Z  transform 
corresponding  to  the  sampled  output  of  a  continuous  filter, 
considering  the  signal  to  be  the  impulse  response  of  the 
filter,  h(t) . 

In  general  there  are  many  known  procedures  for  the  design 
of  continuous  filters  to  meet  various  types  of  specifications. 
These  procedures  are  usually  described  in  the  s  domain  or  in 
the  frequency  domain,  s  =  jco.   The  Z  transform  is  used  to 
translate  continuous  time  realizations  into  discrete  time 
realizations  or  digital  filters  - 

In  control  theory,  when  the  differential  equations  of 
the  systems  are  known,  it  is  often  necessary  to  simulate 
continuous  systems  on  a  digital  computer.   One  approach  is  to 
consider  the  impulse  response  which  is  the  inverse  Laplace 
transform  of  the  system's  transfer  function.   The  Z  transform 
can  be  applied  to  this  signal  after  sampling,  yielding  a  discrete 
simulation  of  the  system. 

Equation  (2.8)  assures  that  the  Z  transform  of  a  digital 
filter  or  system,  when  expanded  in  terms  of  inverse  powers 
of  z,  represents  identically  the  impulse  response  of  the 
continuous  filter  or  system  sampled  at  intervals  T. 

If  X(z)  is  the  Z  transform  of  the  input  signal  to  a  system, 
and  if  H(z)  is  the  Z  transform  of  the  system  transfer  junction, 
then 

Y(z)  =  X(z)  H(z)  (2.23) 
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Here  Y(z)  represents  the  Z  transform  of  the  output  signal, 
and  corresponds  to  the  numerical  convolution  of  the  input 
signal  an  the  system  transfer  function,  as  is  demonstrated 
below 

00 

X(z)  =  Z    x(mT)  z"m 
m=0 

(2.24) 

H(z)  -  Z   h(kT)  z~k 
k=0 

where  h(kT)  is  the  sampled  impulse  response  of"  the  system. 
Substituting  (2.24)  into  (2.23) 

OO       CO 

Y(z)  =  I         Z    x(mT)  h(kT)  z"(k+m)  (2.25) 

m=0  k=0 

and  replacing  subscripts,  such  that  k+m  -  n  or  k  =  n-m,  then 
the  limits  of  the  summations  are  converted  as  follows: 
k  =  0  becomes  n  =  m  and  k  =  °°  becomes  n  =  °°.. 
Thus  (2.25)  becomes 

OO       00 

Y(z)  =  Z    Z    x(mT)  h((n-m)T)  z~n  (2.26) 

m=0  n=m 

Changing  the  order  of  the  summations  and  neglecting  or 
adding  zero  terms  of  the  type  h(-kT)  (k=l,2, .)  yields 

oo       n 

Y(z)  =  Z      Z    X(mT)  h((n-m)T)   z  n  (2.27) 

n=0    m=0 


Z   y(nT)  z  n 
n=0 
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Hence  the  coefficient  y(nT)  of  z    in  Y(Z)  is  given  by 

n 
y(nT)  =  I         x(nT)  h((n-m)T)  (2.28) 

m=0 

Equation  (2.28)  is  the  numerical  convolution  of  the  input 
samples  and  the  digital  filter  sampled  impulse  response. 

The  output  signal  transform,  Y(z),  when  expanded  in  a 
series  of  inverse  powers  of  z  represents  the  output  signal  at 
different  times.   If  X(z)  and  H(z)  are  obtained  from  sampling 
a  continuous  input  signal  and  the  impulse  response  of  a  con- 
tinuous system  or  filter  respectively,  then  the  coefficients 
in  the  expansion  of  Y(z)  will  correspond  to  the  sampled  output 
of  the  continuous  system  or  filter.   Hence  when  synthesizing 
a  digital  filter  using  the  Z  transform  of  a  continuous  filter 
impulse  response  sampled  at  intervals  T,  the  resulting  output 
signal  is  an  exact  description  in  the  time  domain  in  the  sense 
that  the  output  values  of  the  digital  realization  will  be  equal 
to  the  values  of  the  continuous  filter  output  signal  sampled 
at  intervals  T. 

Regarding  the  frequency  spectrum  of  the  filters,  consider 
a  continuous  system  with  transfer  function  H(s) .   If  an  input 
signal  E. (s)  is  applied,  the  output  signal  E  (s)  has  a  Laplace 
transform 

E  (s)  =  H(s)  E. (s)  (2.29) 

The  inverse  Laplace  transform  of  E  (s)  is  the  output 
signal  e  (t)  =  £   {E  (s)},  which  now  is  sampled  at  the  rate  of 
1  /  T  samples  per  second,  yielding  the  sampled  output  signal 
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e  * (t)  as  shown  in  Fig.  2.3a.   The  Fourier  transform  of  this 
signal  is  the  desired  output  spectrum 

E  *(joj)  =  J   {e  *(t)  }  (2.30) 

o  o 

This  spectrum  could  also  have  been  obtained  by  convolving 
the  Fourier  transform  of  the  sampling  function  l  6  [j  (oo  - 


n= 


no)  )  ]  such  that 
s 


E  *(jaj)  =  E  (juj)  *  £     6[j(a)~nu>  )]  (2..31) 

O  O  o 

n=-oo 


=  I     E  (3  [oo-nco  ]  ) 
o  J      s 
n=-°° 


This  signal  processing  may  be  compared  with  the  develop- 
ment of  the  digital  filter  shown  in  Fig.  2.3b.   We  assume  the 
same  sampling  frequency  f   =  1  /  T  and  obtain  the  Z  transform 
H(z)  of  the  filter  function  H(s).   Then  the  Z  transform  of  the 
filter  output  E  (z)  is  given  by 

E  (z)  =  H(z)  E. (z)  (2.32) 

where  E. (z)  is  the  Z  transform  of  the  sampled  filter  input 
e.*(t).   Substituting  z  =  exp(jcoT)  into  E  (z)  yields  the  desired 
output  spectrum  E  *(jco)  .   The  following  question  may  now  be 
asked:   Is  this  output  spectrum  of  the  digital  filter  identical 
with  the  output  spectrum  obtained  in  Eq.  (2.31)?   This  question 
is  discussed  in  the  next  section. 

1.    The  Effect  of  Sampling  on  the  Frequency  Spectrum 

For  a  signal  of  finite  bandwidth  w  radians  called  a 
band  limited  signal  shown  in  Fig.  2.4a  Shannon's  Sampling 
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Theorem  [8],  requires  a  sampling  rate  equal  to  at  least  twice 
the  maximum  frequency  w  of  the  signal  in  order  to  recover  the 
signal  without  distortition  after  passing  it  through  an  ideal 
low  pass  filter.   The  spectrum  of  the  signal  sampled  at  a 
higher  rate  than  the  Nyquist  rate  is  shown  in  Fig.  2.4b..   If 
the  sampling  frequency  is  lower  than  the  limit  set  by  Nyquist, 
then  some  overlapping  of  spectrum  will  occur.   This  is  known 
as  the  aliasing  effect  and  is  illustrated  in  Elg.  2.4c.   It  has 
been  shown  by  L.  J.  Fogel  [9],  that  it  is  possible  to  recover 
a  signal  sampled  at  a  lower  rate  than  Nyquist  (or  one  for  which 
aliasing  has  taken  place)  if  for  every  sample  additional  infor- 
mation about  the  signal  is  known.   For  example,  this  is  possible 
if  its  derivatives  with  respect  to  time  are  known.   However, 
this  idea  cannot  be  applied  here  because  we  restrict  ourselves 
to  single  input  filters  where  only  the  values  of  the  signal  are 
known . 

To  answer  the  previous  question  regarding  comparison 
of  the  two  spectra  of  Fig.  2.3,  consider  the  continuous  case 
of  Fig.  2.3a,  where  for  s  =  jco 

E  (ja>)  =  H(jco)  E.  (ju>)  (2.33) 

After  sampling  this  spectrum,  E  *  ( joo)  is  the  convolu- 
tion of  E  (jco)  and  the  spectrum  of  the  sampling  signal,  which 
is  a  train  of  impulses  separated  by  u   on  the  w  axis  as  given 
by  Eq.  (2.31).   The  effect  is  to  make  the  sampled  spectrum 
repetitive  with  period  w  .   From  Shannon's  Sampling  Theorem, 
the  sampled  output  can  be  reconstructed  without  aliasing  distor- 
tion provided  E  (jco)  is  band  limited  and  the  sampling  frequency 
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(b)   Spectrum  of  a  sampled  band  limited 
signal,  f   >  2w 
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(c)   Spectrum  of  a  sampled  band  limited 
signal,  2jrf   <  2w 


Fig.  2.4.   Effect  of  sampling  frequency  on  the 
frequency  spectrum  of  the  sampled 
signal . 
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is  higher  than  two  times  the  maximum  frequency  of  E  (jw) .   For 

E  (joo)  to  be  band  limited  either  E.  (joj)  or  H(jco)  have  to  be 
o  1 

band  limited. 

Consider  now  the  discrete  case  shown  in  Fig.,  2.3b, 
where  e. (t)  is  sampled  and  the  Z  transform  is  obtained.   From 
Eq .  (2.8) 

(2.34) 


E±(z)  =  E±*(s) 


sT 
z  =  e 


Substituting  Eq .  (2.34)  into  equation  (2.32)  yields 

Eo*(ejwT)  =  H*(ejwT)  E±* (ejuT)  (2.35) 

»  tcoT 

As  may  be  seen,  the  outpt.t  spectrum  E  *  (eJ   )  is  the 

product  of  the  two  spectra  H*  (e    )  and  E.*(e    )  .   Clearly 

both  these  spectra  must  be  undistorted  (no  aliasing)  if" 

E  *  (eJ   )  is  to  be  undistorted.   Hence  for  this  case  the 
o 

requirement  is  different  from  the  case  of  Fig.  2.3a.,   In  view  of 
the  above  remarks  we  conclude  that  it  is  necessary  to  have 
E.  (jco)  and  H(ju))  band  limited  and  furthermore,  one  must  sample 
at  a  frequency  above  the  Nyquist  rate  for  both  spectra  of 
Figs.   2.3a  and  2.3b  to  be  identical. 

As  an  example  consider  the  design  of  a  digital  low 
pass  filter  with  a  cut  off  frequency  of  2  KHz  which  is  required 
to  filter  a  signal  of  10  kHz  bandwidth.   Such  a  filter  must 
be  designed  with  a  sampling  rate  greater  than  20  kHz..  Note 
that  if  the  signal  had  been  obtained  from  a  continuous  filter, 
it  could  be  recovered  by  sampling  at  a  rate  of  only  4  kHz. 

It  is  clear  from  the  above  remarks  and  the  example  given, 
that  there  is  a  maximum  frequency  that  can  be  handled  by  a 
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digital  filter  in  real  time  operation.   In  fact  this  limit  will 
depend  on  the  computing  speed  of  the  processor  and  the  general 
complexity  of  the  filter.   It  is  desirable  to  keep  digital 
filter  designs  as  simple  as  possible,  although  considerations 
of  stability  may  modify  these  considerations. 

2 .    The  Interelationship  Between  the  s  and  z  Planes 

The  interelationship  betv/een  the  s  plane  and  the  z 
plane  can  be  derived  from  the  definition  of  the  variable  z 
namely 

z  =  eST  (2.36) 

A  point  in  the  s  plane,  namely  (a+joo)  will  be 
mapped  into  the  point  in  the  z  plane 

z  =  e<a+jU)T  =  e0T  eJu,T  (237) 

and  so    | z |  =  e 

Therefore,  the  points  in  the  s  plane  with  constant 

damping,  i.e.,  a  =  constant,  map  onto  a  circle  of  constant 

oT 
radius,  e   ,  in  the  z  plane.   Points  in  the  s  plane  with 

constant  frequency,  i.e.   u)  =  constant,  map  onto  radial  lines, 

arg  z  =  wT . 

Furthermore,  for  points  in  the  left  half  s  plane 

a  <  0 ,  and  so 

| z |  =  eaT  <  1   for   T  >  0  (2.38) 

Equation  (2.38)  indicates  that  the  entire  L-H.P. 
maps  into  the  inside  of  the  unit  circle  in  the  z  plane. 
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Although  useful,  the  foregoing  discussion  may  be 
misleading  unless  it  is  applied  properly.   The  misuse  of 
Eq.  (2.36)  occurs  often  in  the  literature.   Take  for  example 
the  paper  by  C.  J.  Greaves  and  J.  A.  Cadzow  [10].   Theijr  mis- 
take consists  of  replacing  the  variable  s  in  the  Laplace 
transform  of  the  continuous  filter  by  its  corresponding  z,  with 
s  =  (1/T)  In  z  as 


H(z)  =  H(s) 


(,2.39a) 
s  =  (1/T)  In  z 


This  equation  is  in  error  and  should  be  replaced  by 


H (z)  =  H*  (s) 


(2.39b) 
s  =  (1/T)  In  z 


The  reason  is  clear  from  the  closed  form  expression 
for  the  Z  transform,  Eq.  (2.18),  where  before  the  substitution 
of  variables  is  made,  a  complex  integral  has  to  be  evaluated. 

However,  there  is  a  direct  connection  between  the 
Laplace  transform  of  the  continuous  function  and  the  Z  trans- 
form of  the  sampled  function.   Namely  the  poles  of  the  Laplace 
transform  of  the  continuous  function  are  mapped  into  the  z  plane 
of  the  sampled  function  using  the  mapping  function  of  (2.36). 

This  can  be  seen  from  equation  (2.22)  where  every 
term  in  the  summation  will  have  a  pole  in  the  z  plane  at 

z  =  e  nT  (2.39c) 

which  indicates  the  mapping  of  the  singularity  X      according  to 
the  mapping  relation  of  (2.36). 
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In  summary  Eq .  (2.36)  may  only  be  used  to  transform 
poles  of  the  continuous  function  in  the  s  plane  to  pole 
locations  of  the  sampled  function  in  the  z  plane.   It  does  not 
represent  a  change  in  variable  except  when  applied  directly 
to  a  sampled  signal.   As  noted  this  subtlety  has  proven  to  be 
a  source  of  many  mistatements  in  the  literature. 
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III.   THE  ALGEBRAIC  SUBSTITUTION  METHOD 

s 

An  equivalent  digital  representation  of  a  continuous 
system  or  filter  is  often  required  in  practice.   This  could 
arise,  for  example,  in  a  digital  simulation  of  a  continuous 
control  problem,  or  in  an  attempt  to  extend  continuous  filter 
theory  to  provide  digital  filter  realizations.   This  latter 
application  is  very  appealing  because  continuous  filter  theory 
is  very  well  developed  and  procedures  exist  to  match  given 
specifications  to  filter  relizations .   Designs  have  been 
tabulated  and  are  available  in  many  handbooks.   If  a  direct 
method  to  obtain  the  equivalent  digital  realization  of  a  contin- 
uous  filter  exists,  it  would  make  all  previous  knowledge  of 
continuous  filter  theory  available  for  the  discrete  case. 

It  is  convenient  when  implementing  a  digital  filter,  to 
obtain  its  representation  as  a  Z  transform,  and  particularly 
as  a  ratio  of  two  polynomials  of  inverse  powers  of  z ..   Such  a 
representation  indicates  the  linear  combinations  which  have  to 
be  performed  upon  past  values  of  input  and   output  signals 
samples,  to  obtain  the  present  value  of  the  output. 

The  most  common  description  of  continuous  filters  in  the  s 
domain  is  as  a  ratio  of  two  polynomials  in  s  and  so  it  is  con- 
venient to  have  a  rational  function  of  z  to  directly  replace 
for  s  in  the  continuous  type  description.   This  is  if 

s  =  f(z)  ((3.1) 

then   H(s)  =  H(f(z))  =  G(z)  .  (3.2) 
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If  the  function  f (z)  is  a  ratio  of  two  polynomials  of 
inverse  powers  of  z  and  if  H(s)  is  a  ratio  of  two  polynomials 
in  s,  then  G(z)  in  Eq.  (3.2)  is  assured  to  be  the  ratio  of  two 
polynomials  in  inverse  powers  of  z. 

This  procedure,  described  in  Eq.  (3.2)  is  called  the 
algebraic  substitution  technique,  and  converts  rational 
polynomials  in  s  into  the  desired  rational  polynomials  in  z. 

It  will  be  shown  that  the  selection  of  a  specific  f (z) 
to  substitute  for  s,  corresponds  to  the  adoption  of  a  numerical 
integration  algorithm  to  integrate  the  phase  state  variables 
of  the  filter  or  system  in  discrete  time. 

A.    MEANING  OF  THE  ALGEBRAIC  SUBSTITUTION  S  =  f (z) 

In  order  to  find  functions  that  can  replace  the  variable 
s  in  the  description  of  the  continuous  filter,  H(s),  the 
meaning  of  the  substitution  must  be  understood.   In  the  first 
place,  the  variable  s  of  the  Laplace  transform  theory  represents 
an  operator  in  the  time  domain,  i.e.,  it  is  a  differentiator. 
Similarly,  the  inverse  of  s,  s   ,  represents  an  integrator. 
Consider  a  filter  whose  transfer  function  is 

a   +  a,  s  +  .  .  .  +  a   s 

H(s)  =  — - ,  where  m  >  n       (3.3) 

b  +  b,  s  +  .  .  .  +  b   s 
o    1  m 

If  this  transfer  function  is  the  ratio  of  the  output  signal 
transform  to  the  input  signal  transform  with  zero  initial 
conditions,  then  Eq.  (3.3)  can  be  expanded  as  follows 

u(5)    -   Y<s)  -  W(s)  .  Y(s)  ,3  4) 

H(s)  ~  xT¥T  "  xTiT   wW  (  } 
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The  following  equations  can  be  formed  arbitrarily 

aisi--  ^ (3.5, 

X(s)     b  +  b,  s  +  b0  s   +  ...  +  b   s 
o     1       2  m 


|{§i=  aQ  +  ax  s  +  ...  +  an  sn  (3.6) 

Rearranging  Eq.  (3.5)  and  (3.6)  yields 

m-1 

b   sm  W(s)  =  X(s)  -  Z    b.  s3    W(s)  (3.7) 

m  .  A   i 

3=0   J 


n 
Y(s)  =  I        a.  sD  w(s)  (3.8) 

j=0  i 

Equations  (3.7)  and  (3.8)  can  be  illustrated  as  in  Fig. 
3.1,  indicating  a  possible  implementation,  as  is  usually  done 
in  an  analog  computer  simulation  for  low  order  systems. 

It  has  been  shown  that  a  transfer  function  of  the  type 
described  in  Eq.  (3.3),  can  always  be  expanded  in  the  form  shown 
in  Fig.  3.1.   From  this  illustration  it  can  be  seen,  that  the 
replacement  of  the  integrations  (denoted  by  1/s  in  the  state 
space  expansion  of  a  continuous  filter)  by  a  function  of  z, 
l/f(z),  corresponds  to  a  direct  algebraic  substitution  of  the 
variable  s  by  the  function  s  =  f (z)  in  the  closed  form  expression 
of  Eq.  (3.3) . 

Furthermore,  the  linear  operation  l/f(z),  should  perform 
an  approximate  operation  in  the  discrete  time  domain  similar 
to  the  operation,  1/s,  of  the  transform  domain,  if  the  responses 
obtained  with  the  digital  filter  are  to  resemble  the  responses 
obtained  from  the  continuous  system. 
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In  view  of  the  above  remarks,  it  can  be  concluded  that 
the  algebraic  substitution  technique  represents  the  adoption 
of  a  numerical  integration  algorithm  to  perform  the  operation 
of  integration  for  the  state  variables  of  the  filter. 

B.    OBTAINING  RECURSIVE  INTEGRATION  ALGORITHMS 

To  derive  recursive  integration  algorithms,  it  must  be 
decided  whether  or  not  the  present  value  of  the  input  is 
available  to  compute  the  present  value  of  the  output.   This 
question  can  be  settled  if  some  assumptions  are  made.   For 
instance,  if  the  time  required  to  perform  the  computations  is 
considered  negligible,  that  is  if  the  time  elapsed  between  the 
present  input  and  the  corresponding  output  of  the  filter  is  of 
no  importance,  then  the  algorithms  can  be  derived  using  inter- 
polation formulas.   On  the  other  hand,  when  this  factor  is 
considered  important  and  the  output,  at  the  present  time  is 
computed  without  the  present  input  value,  then  some  extrapola- 
tion formula  must  be  used.   This  latter  condition  involves  some 
form  of  a  predicting  method,  and  intuitively  seems  to  be  less 
accurate.   Among  several  different  criteria  to  find  recursive 
integration  algorithms,  the  following  methods  are  described 
because  they  lead  to  well  known  integration  formulas  which  are 
used  in  practice. 

a.    The  Adams-Moulton  Method. 

This  method  is  derived  from  the  expression 


*k  =  yk-l  + 


K   f  (t)  dt  (3.9) 

t    X 
rk-l 
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where  f  (t)  is  a  polynomial  of  (n-1) 
x 


-st 


degree  approximating 
the  input  signal,  based  on  the  last  n  available  input  values 
The  recursive  integration  algorithms  for  the  first 
two  degrees  of  approximating  polynomials  are 


Interpolating  Form 
n   (Using  Last  input) 

1  yk=yk-if!  <xk+xx-i> 


Extrapolating  Form 
(Not  using  last  input) 


yk=yk-i+2(3xk-rxk-2 


2   yk=yk-l+r2<5Xk+8Xk-l-Xk-2'    V**-!4!!  <23Xk-l-16Xk-2  +  5Xk-3' 


From  these  algorithms  it  can  be  seen  that  the  first 
degree  approximation  using  the  last  input  value  is  the  well 
known  Trapezoidal  rule  of  integration. 
b.    The  Milne's  Method. 

This  method  correspond  to  the  numerical  integration 
criterion 

-t, 


^k  =  ^k-n  + 


f  (t)  dt 
x 


(3.10) 


'k-n 


-st 
where  f  (t)  is  the  polynomial  of  (n-1)     degree  approximating 

A 

the  input  function  using  the  last  n  available  values  of  the 
inputs . 

Clearly,  when  the  approximating  degree  n  is  unity 
Adams-Moulton  and  Milne's  methods  are  the  same,  giving  identical 
recursive  integration  algorithms.   For  a  polynomial  approxima- 
tion of  second  degree  the  recursive  Milne's  algorithms  are: 
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Interpolating  Form  Extrapolating  Form 

n   (Using  last  input)  (Not  using  last  input) 

2   yk=yk-2+I<Xk+4Xk-l+Xk-2'      ^k=yk-3+£T(2Xk-l-Xk-2+2Xk-3) 


It  can  be  seen  that  Milne's  first  degree  formula  also  leads  to 
the  trapezoidal  rule  and  that  the  second  degree  formula  leads 
to  the  Simpson's  1/3  rule  of  integration,  when  the  last  input 
is  used. 

These  recursive  integration  algorithms  as  well  as  others 
have  been  analyzed  in  the  literature  under  the  general  category 
of  digital  integrators,  as  for  example  in  one  of  the  first 
papers  on  the  subject  which   was  published  by  John  M.  Salzer 
in  1953  [3],  or  in  a  paper  by  A.  P.  Sage  and  R.  W.  Burt  [11]. 

To  obtain  the  frequency  spectra  of  the  algorithms,  one 
must  first  obtain  the  equivalent  Z  transform  representation. 
For  example  consider  the  trapezoidal  rule  of  integration 

yk  =  **-! + 1  (\  +  xk-i»  (3-11) 

The  equivalent  digital  filter  can  be  derived  by  taking 
the  Z  transform 

Y(z)  =  z'1  Y(z)  +  |  (X(z)  +  z"1  X(z))  (3.12) 

and  the  filter  transfer  function  is 

EM  mUsL   =1     i  +  €l  u.i3) 

X  (z)     2    1  -  z  X 
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This  function  corresponds  to  the  algebraic  substitution 

9    1       _1 

s   «   £  -1-  "  z  (3.14) 

T   1  +  Z 

The  frequency  spectrum  of  the  filter  described  in  (3.13) 
is  obtained  by  substituting  z  =  exp(jooT) 


H(eJ   )  =(-)  ^—=  (3.15) 

2   1  -  e  J 


tt  /  JwT)       T        u)T  a-,  1rx 

H(eJ   '  =  -j  ±   cotan  — |  (3.16) 

and  hence 

|H.(ejwT)|  =  cotan  ^  (3.17) 

Arg  H(ejwT)  =  -  j  (3.18) 

These  equations  can  be  compared  with  the  magnitude  and 
phase  of  an  ideal  integrator 

H(s)     =    -  ((3.19) 
s 

H(jtt)    =    tj  (3.20) 

|H(j(o)|   =   i  (3.21) 

Arg   H(joj)    =   -j  (3.22) 

Equation    (3.17)    can   be   expanded  in   a    series   of   the    form 

,„,    jcoT.  .          T      (2           coT         w3T3  ^    ~    1         wT2          r_    _,. 
lH(e         >'    "    2      1ST   "   T"   -ago"-    •••)    -   5"  "12          (3-23) 
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and  therefore  the  magnitude  spectrum  of  the  trapezoidal  rule 
digital  integrator  approximates  the  magnitude  of  an  ideal 
integrator  for  values  of  co  such  that 

2 
-12  <<  w  (3*24) 

or  equivalently 

/12 
w  <<  rf  C3.Z5) 

or 

/3     ~  1 
w  <<   -  oj   ~  ±-  gj  (3.26) 

IT   s    2   s 

In  Salzer's  paper  a  study  of  the  frequency  spectra  of 
several  different  digital  integrators  is  considered  with  a 
comparison  to  the  spectrum  of  an  ideal  integrator. 

Figure  3.2  shows  curves  of  the  trapezoidal  rule  and 
Simpson's  1/3  rd  rule  plotted  against  Q   =  co/oj  .   The  straight 
line  or  real  axis  represents  the  ideal  integrator  and  as  can 
be  seen  the  trapezoidal  rule  falls  below  this  line,  while 
Simpson's  rule  lies  above  this  line. 

In  general  it  can  be  said  that  the  more  complex  the 
integration  algorithm,  the  wider  the  band  of  coincidence  with 
the  ideal  integrator  and  also  the  larger  the  number  of  multi- 
plications to  be  performed  or  the  increase  in  computation 
time.   As  a  result  a  possible  figure  of  merit  might  be  the  ratio 
of  3  db  bandwidth  to  the  number  of  multiplications.   As  a  matter 
of  fact,  simpler  integration  schemes  perform  as  well  or  better 
than  the  complex  ones  in  terms  of  the  above  figure  of  merit. 
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C.    THE  BILINEAR  SUBSTITUTION 

This  particular  case  of  the  algebraic  substitution 
technique  corresponding  to  the  adoption  of  the  trapezoidal  rule 
of  integration  has  a  different  interpretation,  as  is  shown  by 
K.  Steiglitz  [5],   and  also  discussed  by  Gold  and  Rader  [15]. 

In  order  to  avoid  confusion  between  variables,  the 
following  notation  will  be  used: 

s  =  a  +  jw  (3.27) 

(6  +  ja)t  (,3.28) 

£k     —  c 

When  the  bilinear  transformation  is  made,  such  that 

s  =  (§)  ~\  (3.29) 

then  the  entire  jco  axis  is  mapped  onto  the  unit  circle  in  the 
z  plane.   This  can  be  shown  from  the  inverse  relationship  of 
Eq.  (3.29),  namely 

z  =  ±-±i|Lfi  =  1  +  ipg  -  |«|  -  1  C3..30) 

1  -  (2}  S    1  "  (2)jU° 

To  find  the  relationship  between  the  variable  to  in  the 

s  plane  and  the  variable  n     in  the  z  plane  when  s  =  jco,  the 

following  equation  may  be  used 

T 

(—)  -1  out 

z  =  I    ±    K2>     1«  =  eJ-T  =  ej2  tan   =  (i3#31) 

1  -  (f)  j« 
Equation  (3.31)  may  be  reduced  to  the  form 

-n.  =   =■  arc  tan  — —  ((3.32) 
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Equation  (3.32)  indicates  the  transformation  relating  the 
frequency  axis  of  the  s  and  z  planes.   Equation  (3,29)  means 
that  when  the  bilinear  transformation  is  applied  to  a  rational 
expression  in  powers  of  s,  H(s),  a  rational  expression  G(z)  in 
powers  of  z  is  obtained. 

G(z)  =  H(|  f-=rj)  (3.33) 

Furthermore,  when  the  frequency  spectrum  of  the  s  domain 
representation  is  obtained  by  evaluating  the  expression  H(s) 
when  s  =  joo,  and  H(joj)  is  obtained,  and  similarly  if  the  fre- 
quency spectrum  of  the  expression  .in  powers  of  z  is  evaluated 
along  the  unit  circle,  the  function  obtained  is 

G(ej^T)  =  H(|  arc  tan  ^|)  (3.34) 

The  resulting  function  has  the  same  spectrum  as  H(joo)  but 
compressed  into  the  interval  -tt/T  <  -n-  <  tt/T. 

The  nonlinearity  of  the  frequency  scale  can  be  overcome 
by  prewarping  [Ref.  15]  or  frequency  scaling  the  chosen  contin- 
uous realization,  such  that  after  the  transformation  the 
selected  frequency  is  at  the  desired  value.   However,  with  this 
prewarping  technique  only  one  frequency  may  be  scaled  properly 
and  thus  certain  characteristics  of  continuous  filters  as  roll- 
off  per  octave  are  not  preserved.   This  method  of  synthesis  is 
particularly  well  suited  to  synthesize  low  pass  filters  with  a 
cut  off  frequency  much  smaller  than  the  sampling  frequency, 

because  then  the  transformation  function  for  the  frequencies, 

2      -H.T 
cu  =  —  tan  — =■  may  be  considered  linear,  preserving  then  the 

characteristics  of  the  continuous  filter  realization  from  which  the 
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characteristics  were  obtained.   This  situation  can  be  forced 
by  increasing  the  sampling  frequency  or  equivalently  by 
decreasing  the  sampling  interval  T,  and  this  solution  could  be 
the  cure  for  all  problems,  but  this  is  not  so.   As  was 
discussed  in  chapter  two,  reducing  the  sampling  interval,  is 
equivalent  to  reducing  the  time  available  for  computations, 
and  for  a  fixed  computing  speed,  this  results  in  a  reduction 
in  the  maximum  complexity  of  the  filters  that  can  be  implemented. 
Therefore,  the  reduction  in  the  sampling  interval  is  the  least 
desirable  tool  to  be  used. 

Besides  this  consideration  of  computing  time,  there  is 
also  another  consideration  about  the  number  of  bits  required 
in  the  digital  representation  of  the  coefficients  of  the  filter. 

First  it  must  be  noted  that  the  coefficients  in  the  digital 
filters  have  a  finite  accuracy,  and  hence  the  location  of  the 
poles  in  the  z  plane  is  made  on  a  discrete  basis.   To  clarify 
this  consider  a  simple  example  where  the  location  of  a  pole  is 
directly  related  to  the  coefficient  accuracy.   The  filter 

1          z"1 
H(z)  = =  - ^y  (3.35) 

z  -  C-.    1  -  C,  z 

has  a  pole  in  the  z  plane  located  at  the  point  z  =  C-.  .   Because 
of  the  finite  accuracy  in  the  representation  of  C,  this  pole 
can  only  be  located  at  a  finite  number  of  positions.   That  is 
if  a  n  bit  binary  representation  of  the  coefficients  is  used, 
then  the  maximum  accuracy  in  the  pole  position  due  to  quantizing 
effect  is  2 
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Consider  then  the  effect  of  a  reduction  in  the  sampling 
interval  on  the  mapping  of  the  s  plane  into  the  z  plane  under 
the  bilinear  transformation. 

To  visualize  this,  a  grid  in  the  s  plane,  shown  in  Fig. 
3.3,  is  mapped  into  the  z  plane  using  the  bilinear  transforma- 
tion for  values  of  T  =  0.1,  0.05,  0.025,  0.125.   The  results 
are  shown  in  Figs.  3.4,  3.5,  3.6  and  3.7. 

As  can  be  seen,  the  area  in  the  z  plane  into  which,  a 
given  area  of  the  s  plane  is  mapped  decreases  with  T,  the 
sampling  interval.   Conversely  a  given  area  in  the  z  plane 
represents  a  larger  area  in  the  s  plane  as  T  decreases. 

Consider  a  point  z,  in  the  z  plane  distant   one  quantizing 
level  from  the  point  z  =  1  (towards  the  origin  on  the  real 
axis)  and  distant  one  quantizing  level  along  a  vertical  direction 
(above  the  real  axis)  as  shown  in  Fig.  3.8. 

z1  =  (1  -  £n)  +  j  Kn  (3.35) 

The    corresponding   point    in   the    s   plane    is 
2    z.-l  2       -£    +j    £ 

sl  "  t  z~TT  "     f     2-€  +j   5  u.j/j 

1  ^n   J    ^n 


2    5  <-i+ji)  (.2-e  -j  e  ) 

-     — * ^ —  (3.38) 

T  (2-£    )       +    £ 

^n  ^n 


£    2    -,(2.-E    )E  (2-E      -    E)    Cn 

n  n      n      +    j    n    ,      n      „n  (3. .39) 


2  2         J  2  2 

T  (2-£    )       +    C  (2-£    )       +   K 

^n  ^n  sn  n 
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Fig.  3.3.   Grid  on  the  s  plane  to  be  mapped 

into  the  z  plane  using  the  bilinear 
transformation . 
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Fig.  3.4.   Mapping  of  the  grid  of  the  s  plane 
shown  in  Fig.  3.3,  onto  the  z  plane 
using  the  bilinear  trnasformation 
for  T  =  0.1  sec. 
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Fig.  3.5.   Mapping  of  the  grid  of  the  s  plane 
shown  in  Fig.  3.3  onto  the  z  plane 
using  the  bilinear  transformation 
for  T  =  0.05  sec. 
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Fig.  3.6.   Mapping  of  the  grid  of  the  s  plane 
shown  in  Fig.  3.3  onto  the  z  plane 
using  the  bilinear  transformation 
for  T  =  0.025  sec. 
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Fig.  3.7.   Mapping  of  the  grid  of  the  s  plane 
shown  in  Fig.  3.3  onton  the  z  plane 
using  the  bilinear  transformation 
for  T  =  0.0125  sec. 
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Fig.  3.8.   Point  z  =  (  (l-£  ),  j£  )  distant  one 
quantizing  level  from  z  =  I.Q.. 
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2    g  (E  -1)2  E  (1-E  )2 

n   n  ,  .   n     n  -._  .._.. 

=  -   5 o  +  J    p 9  (3.40) 

T    (2-£  )   +  £        (2-C  )   +  5 

n      n         'n     ^n 

b1    ~  ~  -f  +  j  -§   where  ^  <<  1  .  (3.41) 


The  above  derivation  can  be  considered  as  the  mapping  cf 

a  square  of  the  z  plane,  each  side  measuring  one  quantizing 

level,  E  ,  into  an  area  in  the  s  plane  of  side  E   /T.   Because 
n  r  n 

of  quantizing  effects  a  singularity  at  z.,  can  only  be  located 
at  a  finite  number  of  positions,  each  of  which  corresponds  to 
a  position  s . . 

Equation  (3.41)  indicates  that  the  grainniness  of  allowed 
locations  in  the  s  plane  is  dependent  on  the  quantizing  level 
and  the  sampling  interval.   When  using  a  continuous  filter 
design  to  synthesize  a  digital  filter,  the  pole  locations  so 
obtained  may  not  be  conincident  with  these  allowed  positions, 
and  so  it  may  be  necessary  to  shift   them  to  those  allowed 
locations  which  are  closest  to  the  design  locations..   Due  to 
other  considerations,  like  pole  location  sensitivity  in  the  s 
plane,  the  maximum  shifting  allowable  in  the  location  of  a  pole 
in  the  s  plane  may  be  limited  to  a  maximum  distance  A.   If  it 
is  desired  to  make  the  maximum  allowable  shifting  comparable 
to  the  grainniness  in  the  s  plane  the  following  equation  may 
be  established 

A  =   -£  (.3. ,42a) 
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In  the  z -plane  this  grainniness  would  require  n  bits  where 

2"n  =  Sn  (3.42b) 

From  (3.4  2) 

n  =  3.3  log  (~)  (3.43) 

where  n  represents  the  number  of  bits  used  in  the  binary 
representation  of  coefficients  and  A  is  an  indication  of  the 
grainniness  of  the  grid  in  the  s  plane. 

The  lower  bound  of  the  required  number  of  bits  can  be 
determined  from  Eq.  (3.43).   Actually,  the  required  number 
may  be  more,  especially  if  the  poles  are  located  very  far  from 
the  origin  in  the  s  plane,  because  in  this  case  the  distortion 
in  mapping  of  areas  is  bigger  and  can  no  longer  be  neglected 
in  computing  A  through  this  derivation. 

From  the  preceding  discussion,,  it  may  be  concluded  that 
the  finite  representation  of  coefficients  introduces  an  error 
in  the  location  of  the  singularities  of  the  function  under- 
going the  transformation,  and  that  this  error  increases,  that 
is  represents  a  larger  change  in  the  s  plane  when  the  sampling 
interval  is  decreased.   In  other  words,  it  is  not  possible  to 
obtain  a  better  approximation  by  decreasing  only  the  sampling 
interval.   Instead  the  number  of  bits  which  represent  the 
accuracy  of  the  coefficients  must  simultaneously  be  increased 
to  obtain  a  better  approximation  to  the  analog  specification.. 

Because  for  real  time  operation  all  computations  must  be 
performed  within  a  sampling  interval,  the  above  considerations 
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set  a  more  rigorous  limit  than  previously  considered.   This  is 
so  because  an  increase  in  bits  to  represent  the  coefficients 
tends  to  increase  the  computation  time. 

D.    CONCLUSIONS 

The  algebraic  substitution  method  provides  an  easy  method 
to  obtain  a  digital  filter  which  approximates  the  character- 
istics of  a  continuous  filter,  providing  then  the  desired  link 
between  the  continuous  filter  theory  and  the  discrete  filter 
theory,  with  some  restrictions. 

Of  all  the  substitutions,  the  bilinear  transformation  is 
the  easiest  to  apply  and  does  not  increase  the  complexity  of 

4*  Vi 

the  filter,  i.e.,  an    order  continuous  filter  is  transformed 
into  a  n    order  discrete  filter. 

The  simulation  of  continuous  systems  obtained  by  the 
substitution  technique  in  general  can  be  improved  by  decreasing 
the  sampling  interval. 

However,  the  spectrum  obtained  with  a  digital  filter  using 
this  technique  is  only  an  approximation  to  the  spectrum  of  the 
continuous  filter  from  which  it  is  derived,  and  if  the  con- 
tinuous filter  is  itself  an  approximation  to  a  given  specifica- 
tion there  is  no  assurance  that  the  digital  filter  will  also 
be  within  the  specifications. 


51 


IV.   SYNTHESIS  OF  DIGITAL  FILTERS 
FROM  FREQUENCY  SPECTRUM  CHARACTERISTICS 

The  standard  Z  transform  technique  and  the  algebraic 
substitution  method  considered  in  chapters  two  and  three 
represent  an  attempt  to  extend  the  continuous  filter  theory  to 
discrete  or  digital  filter  theory.   There  is  no  method 
available  which  will  guarantee  that  a  digital  filter  will  have 
an  exact  replica  of  the  frequency  spectrum  of  the  continuous 
filter  from  which  it  is  derived,  except  when  the  input  signal 
and  the  filter  transfer  function  are  band  limited  and  sampled 
according  to  Shannon's  theorem  using  the  Z  transform  technique. 
When  using  the  above  methods  to  approximate  the  frequency 
characteristics  of  the  continuous  filters  (which  themselves 
are  an  approximation  to  certain  given  specifications)  there  is 
no  assurance  that  the  resulting  digital  filter  will  meet  design 
specifications  also. 

A  more  direct  approach  to  the  study  of  the  frequency 
spectra  of  digital  filters  seems  to  be  necessary  in  order  to 
provide  better  answers  than  the  techniques  derived  from 
continuous  filter  theory. 

A  part  of  the  material  presented  here  is  an  extension 
of  an  idea  mentioned  by  M.  Martin  [6],   in  which  a  nonrecursive 
filter  is  constructed.   The  approximation  of  Martin  to  a  given 
specification  using  a  Fourier  series  truncated  at  the  n 
harmonic,  requires  2  n  words  of  storage  and  2  n  multiplications 
The  contribution  described  in  this  chapter  leads  to  a  recursive 
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filter,  constructed  from  a  ratio  of  two  finite  Fourier  series 
approximating  a  given  specification.   This  new  procedure 
requires  the  same  storage  and  number  of  multiplications,  but 
allows  one  to  obtain  phase  versus  frequency  characteristics 
other  than  the  linear  phase  versus  frequency  with  slope  obtained 
when  using  Martin's  method. 

The  last  part  of  the  chapter  presents  the  development  of  a 
generalized  technique  to  obtain  recursive  digital  filters 
derived  from  a  magnitude  squared  criterion.   The  magnitude 
squared  function  is  expanded  as  a  ratio  of  two  finite  Iburier 
series.   For  the  approximating  series  to  involve  up  to  the  n 
harmonic,  this  procedure  requires  only  n  words  of  storage, 
which  represents  a  substantial  improvement  over  the  magnitude 
•methods  discussed  in  the  early  part  of  the  chapter. 

A.    SYNTHESIS  OF  NONRECURSIVE  FILTERS  FROM  A  GIVEN"  MAGNITUDE 
SPECTRUM 

In  certain  applications  only  the  magnitude  spectrum  of  the 
output  signal  is  important  and  a  time  delay  of  some  sampling 
intervals  can  be  tolerated.   The  procedure  that  follows  allows 
the  given  magnitude  function  to  be  approximated  with  a  finite 
summation  of  terms  of  a  Fourier  series.   This  procedure  is 
mentioned  by  M.  Martin  [6]   with  reference  to  data  transmission 
filters . 

Consider  the  elementary  filter  whose  transfer  function  is 
of  the  form 

H(z)  =  |  U+z"2]  (4.1) 

A  flow  diagram  of  the  filter  is  shown  in  Fig.  4.1a.. 
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(b)  The  elementary  filter  of  n    order 


Fig.  4.1.   Elementary  filters  used  for  the 

Fourier  series  approximating  technique 
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The  elementary  filter  described  has  a  spectrum  of  the 
form 

H(eJ   )  =  ^[1   +   e  J    ]  (4.2) 

where  z  =  exp(jwT)  has  been  substituted  in  Eq.  (4.1)..   Equation 
(4.2)  can  be  reduced  to  the  form 

H(ejwT)  =  cos  wT  e~ja)T  (4.3) 

and   hence       |H(eJ       )  |    =    |  cos    coT  |  (4.4) 

Arg      H(e^T)    =    -    u>T  (4.5) 

The  elementary  filter  described  in  (4.1)  is  considered 
to  be  of  order  one  and  thus  a  family  of  elementary  filters  of 
higher  order  can  be  generated  in  the  form  of 


Hn(z)  =  |  [1  +  z"2n]  (4.6) 

Such  filters  have  frequency  spectra  of  the  form 

H  (eJ   )  =  cos  nwT  e  J  (4.7) 

n 


|H(ejwT) |  =  |cos  nwT|  (4.8) 

Arg  H(ejwT)  =  -  nooT  (4.9) 

The  filter  of  order  n  is  illustrated  in  Fig.  4.1b. 

Using  these  elementary  filters  it  is  possible  to  approxi- 
mate any  given  magnitude  versus  frequency  characteristic 
|G(eJ   ) | ,    as  will  now  be  shown. 

First  it  must  be  noted  that  if  the  output  of  any  of  these 

filters  described  in  Eq.  (4.6)  is  delayed  m  sampling  intervals, 
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the  frequency  response  becomes  of  the  form 

Hn(z)    =   |    [1    +    z"2n]     '    z"m  (4.10) 


H    (eJ      )    =   cos   nooT  e    J  e    J 

n 

m      -j  (n+m)  00 T  „.  .   ,  n . 

=   cos    nooT   e    J  (4.11) 

Comparing  Eq.  (4.11)  with  Eq .  (4.7),  it  is  clear  that  the 
magnitude  of  the  frequency  spectrum  does  not  change  by  delaying 
the  output.   Only  its  phase  changes  in  the  form  of  an  increase 
in  the  slope  of  the  linear  phase  versus  frequency 
characteristic . 

Therefore,  by  delaying  the  output  signal  of  these  filters 
a  convenient  number  of  sampling  intervals,  it  is  possible  to 
make  their  output  signals  have  the  same  phase  characteristics, 
and  hence  if  several  of  these   specially  delayed  filters  are 
added  in  parallel,  the  resulting  magnitude  spectrum  is  the  sum 
of  the  individual  magnitude  spectra  because  all  the  complex 
quantities  to  be  added  have  the  same  phase. 

An  example  of  the  above  paragraph  is  illustrated  in  Fig. 
4.2a  where  the  elementary  filters  of  zero  and  first  order  are 
delayed  such  that  their  phase  versus  frequency  characteristic 
become  equal  to  the  second  order  elementary  filter..   In  Fig. 
4.2b  the  canonical  form  of  the  filter  is  shown,  indicating 
symmetry  in  the  coefficients. 

Generalizing,  a  filter  synthesized  as  a  summation  of  the 
delayed  elementary  filters  in  the  form 
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X(z) 

o- 


-1 

z 

z 

-1 


-1 


-1 


(a)   Equation  form 


1 


-3_„-l 


l 

2 


HziiK+>-t> 


-j £M+>— 


c,         „    Y(z) 


{> 


Y(z) 


(b)    Canonical    form 


Y(Z)    -    C2    .    Cl      -1  -2         cl      -3         c2      -4 

xTiT  --T+-2Z       +C0Z       +TZ       +TZ 


Fig.  4.2.   Synthesis  of  nonrecursive  filter 

approximating  a  given  spectrum  function 
involving  up  to  the  second  harmonic  of 
the  Fourier  Series. 
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N 
G(z)  =  Z     C   [1  +  z"2n]  '  z"N+n  (4.12) 

has  a  frequency  spectrum  which  can  be  obtained  by  substituting 
=  exp(jcoT)  into  Eq .  (4.12),  yielding 


N 
G(eJb)T)    =      E     C   cos  nwT   e_:,NwT  (4.13) 

n=0    n 

■inT  N 

and   hence         I G (eJ       ) I    =    |Z  C      cos   nwT I  (4.14) 

n=0 
Arg   G(eDa)T)    =    -   NwT  (4.15) 


for  real  coefficients,  C  ,  with  (n=0 , 1, . . . , N) . 

Equation  (4.14)  indicates  the  possibility  of  approximating 
any  given  magnitude  spectrum  with  a  Fourier  series  of  finite 
numbers  of  terms.   It  is  known  that  the  magnitude  spectra  are 
even  functions,  and  repetitive  with  period  oj   =  2tt/T.   This 
characteristic  assures  that  it  will  always  be  possible  to 
expand  a  magnitude  spectrum  function  as  a  Fourier  series  with 
cosine  terms  only,  assuring  also  the  generality  of  Eq.  (4.14). 

Equation  (4.15)  indicates  that  the  resulting  function 

will  be  delayed  N  sampling  intervals,  N  being  equal  to  the 

number  of  harmonics  considered  in  the  Fourier  series  expansion. 

If  F    C   cos  ncoT  >   0   for  all  co  (4.16) 

n=0   n 

then  the  phase  given  by  Eq.  (4.15)  will  not  display  any  dis- 
continuity jumps  and  will  remain  linear.   This  can  always 
be  accomplished  by  proper  adjustement  of  the  coefficient  CQ . 
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It  is  well  known  that  the  Gibb ' s  phenomenon  occurs  near 
a  discontinuity  of  a  function.   In  other  words  no  matter  how 
many  harmonics  are  considered  in  the  Fourier  series  approxima- 
tion, the  summation  of  terms  does  not  converge  to  the  value 
of  the  function  near  a  discontinuity,  but  does  converge  to  an 
overshoot  and  a  undershoot  of  approximately  10%  of  the 
discontinuous  jump.   In  order  to  make  certain  that  (4.16)  is 
satisfied,  the  value  of  C„  must  be  adjusted  to  include  the 
undershoot  of  the  Gibb's  phenomenon.   In  the  case  of  an  idefil 
low  pass  filter  the  Gibb's  phenomenon  is  about  10%  of  the 
discontinuity.   If  Cn  is  adjusted  so  that  the  Fourier  series 
is  always  positive,  this  implies  that  the  attenuation  in  the 
stop  band  is  approximately  2  0  db. 

Some  other  methods  may  be  used  to  obtain  a  greater 
attenuation.   For  instance  the  Lanczos  and  Fejer  methods  to 
attenuate  or  eliminate  the  overshoot  and  ripple  of  the  Fourier 
series  can  be  used  as  discussed  by  R.  W.  Hamming  [12] . 

A  different  approach  to  get  around  this  inconvenience  is 
to  work  numerically,  with  the  given  magnitude  spectrum  as  a 
set  of  data  points,  through  which  a  continuous  function  is  to 
be  fitted.   The  lack  of  discontinuities  assures  that  Gibb's 
phenomenon  will  not  occur. 

A  simple  example  of  nonrecursive  digital  filter  synthesis 
is  now  presented. 

Synthesize  a  nonrecursive  digital  low  pass  filter  approx- 
imating the  magnitude  function  shown  in  Fig.  4.3,  where  o)c=tt/3T 
and  the  order  of  the  filter  is  chosen  such  that  the  approximation 
involves  the  fourth  harmonic  in  the  Fourier  series  expansion. 
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Representing   the   magnitude    function   as    a   Fourier    series 
yields 

•      m  C  4 

\G(eJb}l)  I    =      -2.     +   Z  C      cos   nuT 

1  '  2.  ,         n 

n=l 


fOJ 


where 


C      =   — 
n        to 


cos   ntoT   doo   = —  sin   ncoT 

to      nT 
o  s  -'■o 


to 


— - —      sin   nTco 
nTco  c 


4co  sin   nTco 

c  c 

to  nTco 
s  c 


For   this    example   to   /co      =    1/6    and   to   T   =   tt/3, 


Therefore    the   coefficients   C      become 

n 


r      =   —     sin   n   tt/3 
n        3  n   tt/3 


CQ  =  0.66666   C1  =  0.55133   C2  =  0.27567   C3  =  0.00   C4  =  -0.13783 

Consider  first  the  case  where  C  is  not  adjusted  so  that 
(4.16)  is  not  satisfied. 


jooT 


|G(eD       )|    =    0.33333    +    0.55133    cos    toT   +    0.27567    cos    2toT 


+    0.0    cos    3toT   -    0.13783    cos    4toT 
and    from  Eqs.     (4.14)    and    (4.12) 

m    ^          n    oooo?      "4            0.55133     M            -2,    -3         0.27567 
G(z)    =    0.33333    z         +      ^ [1    +    z      ]z         +   ~ 


M            -4,    -2         0.13783     ri  -2, 

[1    +    z       ]z         -   ~ [1    +    z       ] 
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and  rearranging  terms  the  filter  becomes 

G(z)  -  -  0.068915  +  0.13783  z~2  +  0.27567  z~3  +  0.33333  z"4 

+  0.27567  z"5  +  0.13783  z~6  -  0.068915  z~8 

In  Figs.  4.4  and  4.5  the  magnitude  and  phase  versus 
frequency  characteristics  are  shown  for  the  filter  of  the 
example,  where  the  coefficients  are  obtained  directly  from  the 
Fourier  series  expansion.   The  discontinuities  in  the  phase 
function  can  be  seen  to  correspond  to  every  cross  over  point 
shown  in  the  magnitude  function  by  an  increase  in  attenuation. 
The  lobes  in  the  attenuation  band  of  the  magnitude  function 
are  the  ripples  of  the  Fourier  series  but  distorted   by  the 
logarithmic  scale.   If  C„  is  adjusted  so  that  Cn  =  0.46030  then 

G(z)  =  -  0.068915  +  0.13783  z~2  +  0.27567  z~3  +  0.46030  z~4 

+  0.27567  z~5  +  0.13783  z~6  -  0.068915  z~8 

The  magnitude  and  phase  response  of  this  filter  are  shown 
in  Figs.  4.6  and  4.7. 

In  general  the  foregoing  procedure  requires  2  n  words 
of  storage  for  n  harmonics  considered  in  the  Fourier  series 
approximation.   Furthermore  the  output  signal  has  a  delay  of 
n  sampling  intervals  due  to  the  remarkable  characteristic  of 
this  type  of  filters  of  having  a  linear  phase  with  a  slope 
dependent  only  in  the  number  of  delays  of  the  filter.   Failure 
to  observe  the  condition  stated  in  Eq.  (4.16)  produces   discon- 
tinuities in  the  phase  characteristic  associated  with  each 
crossing  of  the  to  axis  by  the  resulting  Fourier  series.. 
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Fig.  4.4  Magnitude  vs.  normalized  frequency.   Non- 
recursive  filter  with  coefficients  from  Fourier 


series . 
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Fig.  4.5   Phase  vs.  normalized  frequency.   Non- 
recursive  filter  with  coefficients  obtained 
from  Fourier  series. 
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Fig.  4.6    Magnitude  vs.  normalized  frequency.   Non- 
recursive  filter,  coefficients  obtained 
from  modified  Fourier  series  to  obtain 
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Fig.  4.7    Phase  vs.  normalized  frequency.   Nonrecursive 
filter,  coefficients  obtained  from  modified 
Fourier  series  to  obtain  linear  phase. 
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B.    SYNTHESIS  OF  A  SPECIAL  TYPE  OF  RECURSIVE  FILTER  FROM 
A  GIVEN  MAGNITUDE  SPECTRUM  FUNCTION 

The  method  described  in  the  previous  section,  where  non- 
recursive  filters  were  obtained  for  matching  a  given  magnitude 
function,  has  been  generalized  by  the  author  and  alJLows  a 
recursive  filter  to  be  obtained. 

This  special  type  of  filter  has  a  Z  transform  of  the 
impulse  response 


H(z) 


-1              ,            -n    ,                          -2n+l    ,            -2n 
a      +a      ,  z        + .  .  .  +   a  ~  z        + . . . +   a      ,  z                +-  a    z 
n n-1 0 n-1 n 

w       .    v.        -~1  ,    v.      ~m    ,  >    i-  -2m+l    ,    i_      -2m 

b      +b      ,  z        +...+   bz        +...+   b      ,  z  +bz 

ra  m-1  o  m-1  m 

(4.17) 


A  canonical  form  of  this  filter  is  shown  in  Fig.  4.8. 

Note  that  the  polynomials  of  the  numerator  and  the 
denominator  have  an  even  number  of  delays  for  the  input  and 
output  signals,  and  furthermore  are  mirror  image  polynomials 
because  of  the  symmetry  in  the  coefficients.   Such  a  filter 
can  be  considered  to  be  derived  from  the  expression 

H(z)  =  ^r-  (4.18) 


B( 
where  A(z)  = 


-1  -n  -2n+l       -2n 

a   +  a   ,z    +  ...+  a_z    +...+  a   ,z      +  a  z 
n    n-1  0  n-1  n 


(4.19) 


and    B(z)  = 


b   +  b   ,z    +...+  bnz    +...+  b   ,z      +  b_z 
m m-1 0 m-1 m 

(4.20) 
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The  frequency  spectrum  of  such  a  filter  is  given  by 


H(eJuT)  =  Meg, 
B(eJ   ) 


and  from  Eqs.  (4.12) and  (4.13) 


n 
a   +  E    2a,  cos  kwT 

jooT  k=1  e-j(n-m)aiT 

H(eDa)I)  =    (4.22) 

m 
b   +  Z    2b,  cos  kwT 
°    k=l    k 


Equation  (4.22)  shows  the  properties  of  this  special  type 
of  filter  being  considered.   It  indicates  that  its  frequency 
spectrum  is  in  the  form  of  a  ratio  of  two  finite  Fourier  series, 
and  the  phase  is  linear  with  respect  to  frequency. 

When  a  magnitude  versus  frequency  characteristic  is  given 

to  be  approximated  with  a  digital  filter,  if  an  approximation 

in  the  form  of  the  ratio  of  two  finite  Fourier  series  can  be 

obtained,  then  a  method  to  synthesize  the  filter  is  available 

by  using  Eq.  (4.22)  as  will  now  be  explained. 

1 .    The  Expansion  of  a  Magnitude  Function  as  a  Ratio 
of  Two  Finite  Fourier  Series 

Consider  first  the  properties  of  the  Chebyshev 
polynomials,  as  discussed  for  instance  by  R.  W.  Hamming  [13]. 

"The  Chebyshev  polynomials  have  all  the  properties 
of  both  the  Fourier  series  and  the  orthogonal  polynomials;  they 
are,  in  fact,  the  Fourier  functions  cos  n  6  in  the  disguise 
of  a  simple  transformation  of  the  variable 
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If  a  variable  x  is  defined  such  that 

x  =  cos  ojT  (4.23) 

from  where   -1  £  X  _<  1 
then, 

T  (x)  =  cos  (n  arc  cos  x)  =  cos  (ooTri)  (4.24) 

In  view  of   Eq.  (4.24)  an  expansion  in  Chebyshev 
polynomials  in  the  variable  x  corresponds  to  a  Fourier  series 
in  the  variable  toT,  when  the  two  variables  are  related  as  in 
Eq.  (4.23).   This  property  can  be  used  to  approximate  a  given 
magnitude  spectrum  function  as  a  ratio  of  two  finite  Fourier 
series  in  the  following  manner: 

Given  a  function  |G(eJ   ) |   which  has  to  be  approxi- 
mated consider  the  following  transformation 

|G(e^(°T)|  =  G-^ooT)  =  G,(arc  cos  x)  =  H(x)  (4.25) 

If  this  is  done  numerically,  that  is  if  pairs  of 
values  describing  points  of  the  magnitude  versus  frequency 
characteristic  are  given  Eq.  (4.25)  can  be  considered  as  a 
change  in  the  coT  values  of  the  given  points,  maintaining  the 
ordinate  values.   The  interest  in  doing  this  is  that  now  the 
new  function  of  x  can  be  approximated  as  a  ratio  of  two 
polynomials  in  x.   Later  the  polynomials  may  be  expressed  as 
Chebyshev  polynomials  and  then  use  of  the  property  described 
in  (4.24)   leads  to  the  Fourier  series  coefficients.   That  is 
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H(x)  can  be  approximated  by 

P  (x) 
HU)  :  t^~  (4*26) 


and  numerator  and  denominator  polynomials  can  be  represented 
as  an  expansion  of  Chebyshev  polynomials  of  the  form 

n 

E    c   T  (x) 

k=0   K   K 


H(x)  :   m  (4.27) 

E    d   T  (x) 
k=0  K  K 


Using  the  property  of  Eq.  (4.24) 


n 

E  c,  cos  kooT 

•  T     k=0  K 

G(eDwl)|  : (4.28) 

m 

£  d   cos  ktoT 

k=0  K 


and  Eq.  (4.28)  represents  the  approximation  of  a  given  function 
as  the  ratio  of  two  finite  Fourier  series. 

For  synthesis  purposes,  the  coefficients  in  Eq.  (4.28) 
can  be  equated  with  those  of  Eq.  (4.22)  yielding 

ak    =     ck      k  ?    0  (4.29) 

bk    =    dk      k  ^  0  (4.30) 

The  complete  synthesis  procedure  is  illustrated  in 
Fig.  4.9.   The  following  points  should  be  noted: 
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desired  magni 
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Fig.  4.9.   Procedure   to  find  an  approximating  function 
to  a  given  set  of  data  points,  in  the  form  of 
the  ratio  of  two  finite  Fourier  series  with 
cosine  terms  only. 
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a.  If  the  discrete  steps  in  the  abcissae  of 
iGfe-5   )  |  are  uniformly  spaced,  the  abcissae  of  H  (x)  will  be 
nonuniformly  spaced  since  x   =  cos  o>,T. 

b.  Expressing  H(x)  =  P(x)/Q(x)  is  possible  in  many 
ways.   For  the  work  presented  here  IBM  subroutine  ARAT  is  used, 

c.  Expressing  P (x)  and  Q(x)  as  Chebyshev  poly- 
nomials is  performed  using  the  following  tables  taken  from 
McCracken  and  Dorn  [14] . 

1.  Chebyshev  Polynomials 
TQ(x)  =  1 

T1 (x)  =  x 

T2(x)  =  2  x2  -  1 

3 
T. (x)  =  4  x   -  3  x 

T4  (x)  =  8  x4  -7SI  x2  +  1 

Tc  (x)  =  16  x5  -  20  x3  +  5  x 
o 

2.  Powers  of  x  as  functions  of  T  (x) 

n 

1   =  T 
0 

x   =  T. 


X   =  7(T„+TJ 


■2    I,, 

0   2 

3  1 

xJ  =  7T(3T1+T3) 

4  1, 


X   =  7(3T0+4T2+T4) 
XD  =  TT(10T1+5T3+T5) 


The  degree  of  the  numerator  and  denominator 
polynomials  can  be  different,  as  indicated  in  Eq.  (4.30).   The 
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difference  in  polynomial  degree  has  direct  influence  on  the 
phase  characteristics,  as  shown  in  Eq.  (4.22).   Selecting  equal 
degrees  for  numerator  and  denominator,  it  is  possible  to  make 
the  phase  equal  to  zero  for  all  frequencies.   Otherwise  the 
phase  will  change  linearly  with  frequency,  and  the  slope  will 
be  directly  proportional  to  the  difference  in  degree. 

An  example  of  synthesis  using  this  procedure  is  now 
presented.   For  comparison  with  the  M.  Martin  method,   the 
example  given  in  the  previous  section  is  reworked.   The 
magnitude  spectrum  of  Fig.  4.3  will  be  approximated  using  a 
filter  with  eight  delays,  four  in  the  numerator  and  four  in  the 
denominator,  which  implies  the  approximation  of  the  given 
magnitude  spectrum  with  a  ratio  of  two  Fourier  series  involving 
up  to  the  second  harmonic  in  each.   The  method  used  in  the 
example  to  obtain  the  ratio  of  two  Fourier  series  approximating 
the  given  magnitude  spectrum  function,  is  the  numerical  method 
illustrated  in  Fig.  4.9,  using  subroutine  ARAT  of  the  IBM 
scientific  package  to  obtain  the  ratio  of  two  polynomials  in 
the  transformed  variable  x  =  cos  ojT.   The  results  are: 

Given  eight  points  describing  the  magnitude  spectrum 
function  equally  spaced  on  the  cjT  axis  in  the  range 

0  <_  coT  <_   tt  (4.31) 

The  change  to  the  variable  x  is  performed  and  a  ratio  of 
two  polynomials  in  x  approximating  the  points  was  found 

„,  ,    0.01184  +  0.02418  x  +  0.01277  x2  ,.„  „* 

H(x)  =  ~ (4.32) 

1.000  -  1.39494  x  +  0.47740  x 
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The  polynomials  are  now  expanded  in  terms  of 
Chebyshev  polynomials  using  table  shown  previously.   This  yields 

0.01823  T  (x)  +  0.02418  T-  (x)  +  0.00639  T„  (x) 

H(x)  =■-   ii ± i (4.33) 

1.23870  TQ (x)  -  1.39494  T± (x)  +  0.23870  T  (x) 

Using  Eq .  (4.24)  the  ratio  of  Fourier  series  can  be 
obtained. 

ip,  jojT  |    0.01823  +  0.02418  cos  10T  +  0.00639  cos  2coT 
1  [e         '   '    1.23870  -  1.39494  cos  coT  +  0.23870  cos  2wT 

(4.34) 

From  Eqs.  (4.22)  and  (4.28)  it  follows  that 

_,  .    0. 00319+0. 01209z_1+0.01823z~2+0.01209z~3+0.00319z~4 

G(z)  =  — ^^ T3 =4 

0. 11935-0. 69747z   +1.23870z   -0.69747z   +0.11935z 

(4.35) 

The  magnitude  and  phase  of  the  spectrum  has  been 
plotted  and  illustrated  in  Figs.  4.10a  and  4.10b.   As  a  comparison 
with  M.  Martin's  method  shown  in  Fig.  4.4,  it  can  be  noted  that 
there  is  an  increase  in  attenuation  and  elimination  of  lobes 
as  well  as  zero  phase  shift. 

C.    SYNTHESIS  OF  RECURVSIVE  FILTERS  FROM  MAGNTIDUE  SQUARED 
FREQUENCY  SPECTRUM  FUNCTIONS 

In  the  previous  section  the  filter  is  restricted  to  a 

minor  image  coefficient  realization  starting  with  the  magnitude 

of  a  transfer  function.   In  this  section  a  new  approach  is 

developed  to  synthesize  recursive  digital  filters,  from  a  given 

magnitude  squared  function. 
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Fig.  4.10a  Magnitude  vs.  normalized  frequency. 

Recursive  filter,  coefficients  obtained 
from  numerical  approximation  described  in 
Fig.  4.8. 
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Fig.  4.10b  Phase  vs.  normalized  frequency.   Recursive 

filter,  coefficients  obtained  from  numerical 
approximation  described  in  Fig.  4.8. 
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It  is  shown  that  a  digital  filter,  described  in  its  most 
general  form  has  a  magnitude  squared  frequency  spectrum  which 
is  a  ratio  of  two  finite  Fourier  series  with  cosine  terms  only. 
A  method  for  obtaining  an  approximation  to  such  a  function, 
has  been  described  in  the  previous  section,  using  a  change  of 
variable  and  expansion  of  the  approximating  function  in  a 
series  of  Chebyshev  polynomials.   In  this  section  the  coeffi- 
cients of  the  filter  are  obtained  by  equating  the  coefficients 
of  the  Fourier  series  expansion  with  the  coefficients  of  the 
filter  according  to  some  nonlinear  relationships  which  are 
developed. 

1.    Developement  of  the  Theory 

Consider  a  generalized  description  of  a  digital 

filter,  of  the  form  , 

-1        -2  -k   £    a.z"3 

a,,   +   an  z    +  a~  z    +...+  a,  z       ,,   n 

H(z)  =  -2 1 1 * =  lz°_J (4.36) 

-1         _o  _    m 

b„  +  b,  z    +  b0  z    +.  .  .+  b  z  m  I        b.  z"1 
0     12  m      .  0   k 

The  spectrum  of  this  filter  can  be  obtained  by 
evaluating  the  magnitude  and  phase  of  Eq.  (4.17)  for  values  of 
z  along  the  unit  circle 

z  =  e^  =  cos  wT  +  j  sin  wT  (4.37) 

Substituting  (4.37)  into  (4.36),  yields 

-jtoT                   -jcoT    ,                         -jkwT 
.    m           a„    +   a,    e    J         +   a0    e    J         +...+    a,     e    J 
H/pDwT      _      0  1 2 k ,4    38) 

hn   +  bn    e    J         +  bn   e    J         +  ...+   b      e    J 
0  1  2  m 
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Consider 


,  jwT   _  N(jfa)) 
H(S    }  ~  D(ju>) 


|d(jo))  r 


H(eJwT}  |2  =  Re2{N(jo))  }  +  Im2{N(jgj)  }  (4.39) 

Re2{D(jco)  }  +  Im2{D(jaj)  } 


(a  +Z   cos  nwT)   +  (I    a   sin  ntoT) 

|H(e^T)|2  =  _lS=J n=1  " (4.40) 

2     m  2 

(b  +  Z   cos  nwT)   +  (Z    b   sin  nwT) 

°  h=l  n=l   n 

which  can  be  reduced  to  the  form 

k     2   k-1  k-2 

Z   a   +  2Z   a  a  ,,  cos.  wT+2Z   a  a  ,~  cos  2wT  +  .. 
•  mo      ft  n      -.  n  n+i  „  n  n+2 

|H(e3wT)  |2  =  ^ 2=° 5=2 

1        '     m    p   ra-1  m-2 

Z   b   +  2Z   b  b  ,,  cos  u)T+2Z   b  b  ,_  cos  2wT  +.. 

ft  n      ft  n  n+1  n   n  n+2 

n=0      n=0  n=0 

(4.41) 

k 

Z    c   cos  nwT 

|H   J03T  |2  =  n=0_^ 

i        i     m 

Z    d   cos  na3T 
n=0 


Proof  of  Eq.  (4.41)  is  given  in  Appendix  A.   Equation 

(4.41)  indicates  that  the  magnitude  squared  frequency  spectrum 

is  the  ratio  of  two  finite  Fourier  series  with  cosine  terms 

only.   Equation  (4.42)  shows  the  relationships  between  the 

coefficients  of  the  Rmrier  series,  c   and  d   and  the  coeffi- 

n      n 

cients  of  the  filter  a.  and  b..  Thus,  from  (4.22)  and  (4.23) 
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k-n 

c       = 
n 

2 

j=0 

m-n 

a  . 

aj+n 

d     = 

n 

2 

b. 
3 

Vn 

(4.43a) 


(4.43b) 


follows 


This  relationship  can  be  written  in  matrix  form  as 


a0   al    a2 


0    art    a, 
0     1 


0    0 


0    0 


lk-l 


lk-2 


0 


ac 

al 

a^ 

z 

ak 

_ 


1/2 


'2/2 


^Ck/2_ 


(4.44) 


Similar  relationships  hold  for  the  denominator 

coefficients.   As  discussed  in  the  previous  section,  it  is 

possible  to  find  an  approximation  to  an  even  function  as  a 

rational  expression  of  Fourier  series  so  that  c   and  d   are 
^  n      n 

known.   The  coefficients  of  the  Fourier  series,  and  hence  the 
vector  of  c's  in  Eq.  (4.44),  will  be  obtained  from  the  expansion 
of  the  desired  spectrum  function.   Solving  for  the  a's  in 
Eq.  (4.44)  yields  the  coefficients  of  the  numerator  of  the 
digital  filter  realization.   An  identical  procedure  can  be 
followed  for  the  coefficients  of  the  denominator. 

This  set  of  nonlinear  equations  has  several  solutions, 
in  fact,  as  k  gets  larger,  the  number  of  solutions  also 
increases.   This  is  not  surprising  because  of  the  nature  of  the 
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problem.   That  is,  given  a  filter,  a  unique  magnitude  squared 
frequency  spectrum  function  can  be  found,  but  the  inverse 
procedure  is  not  unique.   Given  the  magnitude  of  the  spectrum, 
the  description  of  the  filter  is  incomplete  because  there  is 
no  information  about  the  phase  characteristics. 

Consider  an  example  of  the  use  of  this  synthesis 
method  to  approximate  a  given  magnitude  squared  function.   For 
comparison  with  the  example  given  the  same  magnitude  spectrum 
shown  in  Fig.  4.3  will  be  squared.   What  results  in  the  same 
function,  which  is  now  approximated  with  a  filter  using  only 
two  delays  for  the  output  and  input  signals.   The  nonlinear 
equations  to  be  solved  are  of  the  form 

a0   +  al   +  a2   =  C0  (4.45) 


aQ  ax  +  ax  a2  =  c±/2  (4.46) 


a   a0         =  c0  /0  (4.47) 

S  Z  z/  z 


Solving  this  equations  algebraically  yields 


an  =     /  c   +  c     V  (c   +  c0)2   -   Cl2  C4..48) 

1    —   /   o     2—      o     2       1 


a2  ~   2a 


c2  (4.49) 


C9 

±_  (4.50) 
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Equation  (4.48)  indicates  there  are  four  solutions 
for  a.,  and  Eq.  (4.49)  indicates  there  are  two  solutions   of  an 
for  each  a, .   Hence  there  are  eight  sets  of  coefficients 
satisfying  the  given  nonlinear  equations.   The  same  is  true 
for  the  denominator  set  of  coefficients  and  therefore  there 
are  64  solutions  for  the  filter.   In  general  not  all  the 
solution  may  be  valid  if  only  real  coefficients  are  being 
considered. 

Using  the  results  obtained  in  previous  example, 
Eq.  (4.34),  the  magnitude  squared  function  is  expanded  in  the 
same  ratio  of  Fourier  series,  yielding 

.  ,  jcoT  ,2    0.01823  +  0.02418  cos  a)T  +  0.00639  cos  2coT  ...    . 
1  {e         ;i     1.23870  -  1.39494  cos  ojT  +  0.23870  cos  2coT  (q'^1) 

Substituting  these  coefficients  into  (4.48),  (4.49), 
(4.50),  two  sets  of  solutions  are  obtained,  called  filter  A 
and  B  respectively. 


FILTER  A 

FILTER  B 

ao 

= 

0, 

.038983093 

0, 

.081958605 

al 

= 

0 

.099965521 

0, 

.099965521 

a2 

= 

0 

.081958605 

0, 

.038983093 

b0 

= 

-0 

.284898296 

-0 

.418921424 

bl 

= 

0 

.990978207 

0 

.990978207 

b„ 

_ 

-0 

.418921424 

-0 

.284898296 

The  magnitude  and  phase  versus  frequency  characteristics 
are  shown  in  Fig.  4.11,  4.12,  4.13,  4.14. 
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Fig.  4.11    Magnitude  vs.  normalized  frequency. 

Recursive  filter,  coefficients  obtained 
from  magnitdue  squared  criterion.  Filter 
A. 
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Fig.  4.12    Phase  vs.  normalized  frequency.   Recursive 
filter,  coefficients  obtained  from 
magnitude  squared  criterion.   Filter  A. 
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Magnitude  vs.  normalized  frequency. 
Recursive  filter,  coefficients  obtained 
from  magnitude  squared  criterion.  Filter 
B. 
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Fig.  4.14    Phase  vs.  normalized  frequency.   Recursive 
filter,  coefficients  obtained  from 
magnitude  squared  criterion.   Filter  B. 
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Comparison  of  filters  A  and  B  indicates  the  same 
magnitude  for  both  filters  and  a  different  phase  characteristic 
as  was  expected.   With  these  sets  of  coefficients,  still  two 
more  combinations  are  possible,  that  is  the  numerator  of  A 
and  the  denominator  of  B  and  vice  versa.   At  this  point, 
nothing  has  been  said  about  the  phase  characteristic  of  this 
synthesis  method.   It  is  left  as  a  suggestion  for  further 
research,  the  introduction  of  phase  characteristics  restrictions 
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V.   SUMMARY  AND  SUGGESTIONS 
FOR  FURTHER  RESEARCH 

The  two  main  new  ideas  introduced  in  the  thesis  are: 

1.  A  generalization  of  Martin's  synthesis  technique  to 
include  recursive  filters,  which  allows  filters  with  zero 
phase  or  linear  phase  characteristics,  to  be  obtained. 

2.  A  new  synthesis  technique,  which  depends  upon  the 
expansion  of  a  given  function  as  a  ratio  of  two  Fourier  series. 
The  resulting  technique  involves  the  solution  of  simultaneous 
nonlinear  algebraic  equations,  but  the  filter  requires  half 
the  storage  of  Martin's  method  for  the  same  degree  of  approxi- 
mation.  The  number  of  multiplications  to  be  performed  is  also 
decreased  by  one  half. 

The  foregoing  studies  have  given  rise  to  several  basic 
questions  for  further  research: 

1.    Explore  the  concept  of  expanding  a  given  function 
as  a  ratio  of  two  finite  Fourier  series  directly  from  Fourier 
series  expansion  theory,  as  versus  the  change  of  variable 
technique,  involving  Chebyshev  polynomials,   developed  in  this 
thesis.   It  would  appear  that  an  unlimited  number  of  ratios 
are  possible.   If  this  is  the  case,  filter  stability  can  be 
assured  by  selecting  the  poles  of  the  digital  transfer  function 
to  lie  inside  the  unit  circle  of  the  z  plane,  and  then  selecting 
the  coefficients  of  the  numerator  to  obtain  the  desired 
magnitude  spectrum. 
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2.  Consider  synthesis  procedures  starting  from  a  class- 
ical specification  for  the  filter  in  terms  of  tolerance  bands 
and  attenuation  rates.   Determine  how  these  design  criteria 
can  be  met. 

3.  Find  a  general  method  to  solve  the  nonlinear  equa- 
tions, Eq .  (4.44),  which  were  developed  for  the  filter  synthesis 
proceedure  starting  from  a  magnitude  squared  function  as 
discussed. 


89 


Appendix  A 
Derivation  of  Equation  (4.22) 

k  k 

|H(eDwT)|2      =    (aQ    +    E         aR   cos    nwT)"+(Z         an    sin   nooT)  2  (A.l) 
n=l n=l 

■K.  ~y  JC  ry 

b.    +    Z        b      cos    nooT)~+(Z         a      sin  mooT) 
0  ,       n  ,      m 

n=l  n=l 


consider  the  numerator  only 


N(wT)  I   =  an    +    Z    a   cos  ncoT)   +  (Z    a   sin  nwT) 
1      0     -,   n  ,   n 

n=l  n=l 

=   a~    +2a~Z      a      cos   ncoT+(Z      a      cos    nooT)    +  (Z      a      sin   nooT) 
0  0      •,    n  ,    n  in 

n=l  n=l  n=l 

(A. 2) 
The    last   two    terms    can   be   expanded   as    follows: 


j\.  <-j  JS.  ry  T\.  *y  ry 

(Z      a      cos  nooT)       +    (Z      a  sin   ncoT)       =    Z      a        cos      nwT 

,    n  ,    n  ,    n 

n=l  n=l  n=l 


k-1  k  '  k  2  2 

+    2Z      a      cos    nooT    Z  a.    cos    jwT   +    Z      a        sin      nwT 

n=l   n  j=n+l    3  n=l   n 

k-1  k 

+    2Z      a      sin   nooT    Z  a.    sin    jwT  (A.  3) 

n=l   n  j=n+l    J 

k  2  2  2  k_1        r  k 

Z      a       (cos      nwT+sin      mooT)+2Z      a     jcos    nwT    Z  a.    cos    jwT 

n=l  n  n=l   n  l  j=n+l  3 

k 
+    sin   nooT    Z  a.    sin    jwT^  (A. 4) 


n=n+l 
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k         2  k"1  rk 


E      a        +    2    E        <  E  a      cos    nwT   a  .    cos    ja)T 

n=l  n  n=l  ^  i=n+l   n  3 


k 
+    E  a      cos    ncoT   a.    sin   jtoT^  (A. 5) 

n=n+l    n  ^  J 


k  2  k-1      k 

E      a         +    2    E        {  E  aa.    {cos    ncoT   cos    jcoT 

n=l   n  n=l  ^  j-n+1   n   J 


+   sin   nwT    sin    juiT}  \  (A. 6) 

Use    the    trigonometric    identity 

{cos    nojT   cos    jwT   +    sin   nooT    sin      jooT}    =   cos     (j-n)ooT  (A.  7) 


k  ~  k-1      k 


Ea        +2E       Ie  aa.    cos    (j-n)coT) 

i=l   n  n=l  1  j=n+l   n    ^  J 


changing  subscripts,  i  =  j-n  ,   j  =  n+i   the  limits  of  the 
second  summation  become 

j  =  n+1   ■*■    i  =  1 
j  =  k     -*■     i  =  k-n 


k     „      k-1  k-n 

Ea  +2E    E    aa.  cos  i  coT  (A.  8) 

,  n         ,  .  ,   n  n+i 
n=l        n=l  i=l 


Changing  summation  order  Eq.  (A. 8)  becomes 

k     7  k-1   K-i 

E   a   +  2  I   (  E    a  a   .  I  cos  i  a»T  (A.  9) 

1  n  .  -,  1   t   n  n+i J 

n=l  i=l  k n=l        J 
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Combining  (A. 9)  with  (A. 2)  yields 

9       9  9 

|N(u)T)  I   =  an      +  2  E    aAa.  ,  .  cos  i  uT  +  I    a 
0       ±=1   0  O+i  n=1   n 

k-1   k-1 

+  2  I     £   a  a  , .   cos  i  ojT  (A. 10) 

.  ,     ,  n  n+i 
i=l   n=l 

k     2  k-1 

=  E   a    +2  a_a.  cos  k  wT  +  2  E    a„a0 ,   cos  i  wT 
~  n        Ok  .  ,   0  O+i 

n=0  1=1 

k-1      k-1 

+    2    E  E      a    a       .    cos    i    ojT  (A.  11) 

.    ,  ,    n   n+i 

l-l      n=l 

k  2  k-1   k-i 

=    E      a        +    2    ana,     cos   k   ooT   +    2    E         E      aa,.    cos    i    wT 
r,    n  Ok  .    ,  ~    n   n+i 

n=0  i=l   n=0 

(A. 12) 

k  2  k        k-1 

=Ea        +2E         E         aa,.    cos    i    wT  (A. 13) 

A   n  .    ,         -      n  n+i 

n=0  i=l   n=0 


The  same  relationship  holds  for  the  denominator,  and  so 
established  the  validity  of  Eq.  (4.22). 
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